install.packages("pacman")
pacman::p_load(tidyverse, survival, readxl, ggplot2, survminer)

data <- read_excel("data.xlsx")

# Remove NA values
data$export_city <- as.numeric(data$export_city)
data$disaster <- as.numeric(data$disaster)
data <- na.omit(data[, c("city","population", "democratic", "density", "export_city", "disaster", 
                        "risk", "readiness", "coastal", "t2", "gcom", "airport")])

# Suppressing unnecessary messages
options(warn=-1)

# Declare survival dataset
data$surv <- ifelse(data$gcom == 1 & data$t2 < 11, 1, 0)

# descriptive statistics
summary(data)
summary(data$export_city) / 1000000

data$population <- log(data$population + 1)
data$density <- log(data$density + 1)
data$export_city <- log(data$export_city + 1)

##Analysis

# Fit the Cox model
cox_model1 <- coxph(Surv(t2, surv) ~ readiness, data = data)
cox_model2 <- coxph(Surv(t2, surv) ~ readiness + population  + density + democratic + 
                      coastal + risk + export_city, data = data)
cox_model3 <- coxph(Surv(t2, surv) ~ readiness + population  + density + democratic + 
                      coastal + risk + export_city + frailty(city), data = data)

# PH assumption
PH <- cox.zph(cox_model2, transform="km")
PH

# Extract the summary of the Cox model
cox_summary1 <- summary(cox_model1)
cox_summary2 <- summary(cox_model2)
cox_summary3 <- summary(cox_model3)

cox_summary1
cox_summary2
cox_summary3

extractAIC(cox_model1)
extractAIC(cox_model2)
extractAIC(cox_model3)

## counterfactual scenarios
library(MASS); library(simcf); library(tile)
# for simcf and tile, please refer to: https://faculty.washington.edu/cadolph/index.php?page=60

simCoxPH <- function(sims, before, after, baseline=0) {
  crossRR <- sims^(after - before)
  res <- c(mean(crossRR), quantile(crossRR, probs=c(0.025, 0.975)))
  res
} 

sims <- 10000
simbetas <- exp(mvrnorm(sims, coef(cox_model2), vcov(cox_model2)))

CFnames <- CFs <- NULL

CFnames <- c(CFnames, "Scenario 1")
CFs <- rbind(CFs,
             simCoxPH(simbetas[,1],
                      mean(data$readiness),
                      mean(data$readiness) + sd(data$readiness)
             )
)

CFnames <- c(CFnames, "Scenario 2")
CFs <- rbind(CFs,
             simCoxPH(simbetas[,2],
                      mean(data$population),
                      mean(data$population) + 1
             )
)


CFnames <- c(CFnames, "Scenario 3")
CFs <- rbind(CFs,
             simCoxPH(simbetas[,4],
                      0,
                      1
             )
)

CFnames <- c(CFnames, "Scenario 4")
CFs <- rbind(CFs,
             simCoxPH(simbetas[,7],
                      mean(data$export_city),
                      mean(data$export_city) + 1
             )
)

trace <- ropeladder(x=CFs[,1],
                    lower=CFs[,2],
                    upper=CFs[,3],
                    labels=CFnames,
                    entryheight=0.2,
                    col=c("#e6ab2b", "#8ee62b", "#2b83e6", "#aa2be6"),
                    size=0.5, lex=1.5,
                    lineend="square",
                    plot=1
)

vertmark <- linesTile(x=c(1,1), y=c(0,1), layer=12, plot=1)

file <- "Figure1"

tile(trace, vertmark,
     output=list(outfile=file, width=7.0),
     limits=c(0, 20),
     gridlines=list(type="xt"),
     topaxis=list(add=T, log=F, at=c(1, 5, 10, 15, 20), 
                  labels=c("1x", "5x", "10x", "15x", "20x")),
     xaxis=list(log=F, at=c(1, 5, 10, 15, 20), 
                labels=c("1x", "5x", "10x", "15x", "20x")),
     width=list(plot=2),
     height=list(topaxistitle=2)
)

tile(trace, vertmark,
     limits=c(0, 20),
     gridlines=list(type="xt"),
     topaxis=list(add=TRUE, log=TRUE, at=c(0, 1, 2, 5, 10, 20), 
                  labels=c("0x", "1x", "10x", "15x", "20x")),
     xaxis=list(log=TRUE, at=c(0, 5, 10, 15, 20),
                labels=c("0x", "5x", "10x", "15x", "20x")),
     topaxistitle=list(labels="Relative risk of joining GCoM", y=0.8),
     width=list(plot=2),
     height=list(topaxistitle=2),
     output=list(outfile=file,width=7.0)
) # change the font later

data %>% 
  select(city, readiness) %>% unique() %>% 
  arrange(sort(readiness)) %>% View()

summary(data)